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Abstract 

This paper presents an analysis of the effects of noise and precision on a simplified model of the clarinet driven by a 
variable control parameter. 

In a previous work [5], the dynamic oscillation threshold of a simplified model of a clarinet is obtained analytically. In 
the present article, the sensitivity of the dynamic threshold on precision is analyzed as a stochastic variable introduced 
in the model. A new theoretical expression is given for the dynamic thresholds in presence of the stochastic variable, 
providing a fair prediction of the thresholds found in finite-precision simulations. The dynamic thresholds are found to 
depend on increase rate and are independent on the initial value of the parameter, both in simulations and in theory. 

Keywords: Musical acoustics, Clarinet-like instruments. Iterated maps, Dynamic Bifurcation, Bifurcation delay, Tran- 
sient processes, Noise, Finite precision. 



1. Introduction 

A previous article by the authors [5] analysed the behavior 
of a simplified model of a clarinet when one of its control 
parameters (the blowing pressure) increases linearly with 
time. The oscillations corresponding to the production 
of sound start at a much higher threshold of the blowing 
pressure than the one obtained in a static parameter case. 
The dynamic threshold was described by an analytical 
expression, predicting that it does not depend on the 
increase rate of the blowing pressure (within the limits of 
the theory, i.e. slow enough increases), but that it is very 
sensitive to the starting value of the linear increase. 

These results are reproduced by simulations of the 
model, but only when very high precisions are used in the 
simulations (Fig. 10 in [5]). Running the simulation with 
the normal double-precision of a CPU results in much 
lower thresholds, although higher than the static ones. In 
contrast with the theory and simulations using high pre- 
cision, the threshold depends on the rate of increase of 
the parameter, but doesn't depend on the starting value 
of the blowing pressure. 

The aim of the present article is to estimate the dy- 
namic bifurcation threshold in simulations performed 
with finite precision. The effect of finite numerical pre- 
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cision in simulations is modeled as an ad hoc additive 
white noise with uniform distribution. This hypothesis 
is tested in section 3. In section 4, a mathematical rela- 
tion is derived for the behavior of the model affected with 
noise, following a general method given by Baesens [3]. 
The resulting theoretical expression of the dynamic os- 
cillation threshold is compared to numerical simulations 
and its range of validity is discussed. The clarinet model 
and major results from [5] are first briefly recalled in sec- 
tion 2. 

A table of notation is provided in table I . 

2. Dynamic oscillation threshold of the 
clarinet model without noise 

2.1 Clarinet model 

The instrument is divided into two functional elements: 
the exciter and the resonator. The exciter of the clarinet 
is the reed-mouthpiece system described by a nonlinear 
characteristics relating the instantaneous values of the 
flow u{t) across the reed entrance to the pressure dif- 
ference Ap(f) - PmW - pit) between the mouth of the 
musician and the clarinet mouthpiece [10, 9]. The reed is 
simplified into an ideal spring without damping or iner- 
tia. The resonator is approximated to a straight cylinder, 
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Figure 1: Graphical representation ofj']"'" for different precisions (prec. = 7, 15, 30, 100, 500 and 5000) with respect to the slope e 



th 



and for 70 = 0. Results are also compared to analytical static and dynamic thresholds: jst and jj 



Table 1: Table of notation. All quantities are dimensionless. 



Table of Notation 


G 


iterative function 


r 


musician mouth pressure (control parameter) 


( 


control parameter related to the opening of the 
embouchure at rest 


p+ 


outgoing wave 


p 


incommg wave 


P-* 


non-oscillating static regime of p^ (fixed points of 
the function G) 


<P 


invariant curve 


w 


difference between p^ and (p 


e 


increase rate of the parameter 7 


a 


level of the white noise 


Jst 


static oscillation threshold 


Tdt 


dynamic oscillation threshold 


'dt 


theoretical estimation of the dynamic oscillation 
threshold of the clarinet model without noise or in 
"deterministic" situation 


~th 
'dt 


theoretical estimation of the dynamic oscillation 
threshold of the noisy clarinet model in " sweep - 
dominant" situation 


vth 
'dt 


general theoretical estimation of the dynamic oscil- 
lation threshold of the noisy clarinet model, both 
for a "sweep-dominant regime" or a "deterministic 
regime" 


'dt 


dynamic oscillation threshold calculated on nu- 
merical simulations 



described by its reflection function r(t). Assuming per- 
fect reflexions at the open of the resonator (no radiation 
losses) and ignoring viscous and thermal losses the reflex- 
ion function becomes a simple delay with sign inversion. 

The solutions p{t) and u{f) of the model depend on 
two control parameters: 7 representing the blowing pres- 



sure and f related to the opening of the embouchure at 
rest. In this work, the control parameter ( is always con- 
stant and equal to 0.5. Using the variables p^ - ^[p + u] 
and p~ - i(p-u) (outgoing and incoming pressure 
waves respectively) instead of the variables p and u, the 
state of the system can be calculated at regular intervals 
T - 2Llc, the round trip time of the sound wave along the 
resonator. With these assumptions, the nonlinear system 
becomes an iterated map [12, 13, 11] : 



Pt^G[pl_^,j]. 



(1) 



The function G can be written explicitly for C < 1 (see 
Taillard et al. [15]). When the control parameter 7 is con- 
stant, for low values of 7 the solution of eq. (1) stabilises 
at an equilibrium point, and for a critical value jst a flip 
bifurcation occurs leading to a periodic regime that cor- 
responds to sound production. 

2.2 Dynamic bifurcation 

For a linearly increasing control parameter 7, eq. (1) is 
replaced by eq. (2a) and (2b) : 



Jn^Jn-l + e- 



(2a) 
(2b) 



The theory derived in section 4 requires that the param- 
eter 7 increase slowly, hence e is considered arbitrarily 
small (e« 1). 

Because of the time variation in the control parame- 
ter 7, the bifurcation point corresponding to the appear- 
ance of the oscillations is shifted from the static oscilla- 
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tion threshold jst [6] to the dynamic oscillation threshold 
jdt [5]- The difference Ydt~Jst is called the bifurcation 
delay [3, 7]. The previous article by the authors [5] pro- 
vides an analytical study of the dynamic flip bifurcation of 
the clarinet model (i.e. system (2)) based on applications 
of dynamic bifurcation theory proposed by Baesens [3] . 
The main focus of this work is on the properties of the 
dynamic oscillation threshold, recalled hereafter. 

The trajectory of the system in the phase space (here 
constituted of a single variable p^) through time is called 
the orbit. The dynamic oscillation threshold is defined as 
the value of j for which the orbit escapes from the neigh- 
borhood of the invariant curve (f)(j,e]. This definition 
is different from the one used in [5] where the dynamic 
threshold was defined as the value of gamma for which 
the orbit starts to oscUlate. 

The invariant curve is the nonoscillating solution of 
the system (2) . It plays the role of an attractor for variable 
parameters similarly to the role of the fixed point in a 
static case. The invariant curve is written as a function of 
the parameter, invariant under the mapping (2) and thus 
satisfying the following equation: 



(/)(7,e) = G(0(7-e,e),7). 



(3) 



The procedure to obtain the theoretical estimation y^^ 
of the dynamic oscillation threshold is as follows: a the- 
oretical expression of the invariant curve is found for a 
particular (small) value of the increase rate e (i.e. e « 1). 
The system (2) is then expanded into a first-order Taylor 
series around the invariant curve and the resulting linear 
system is solved analytically. Finally, y^'j is derived from 
the analytic expression of the orbit. 

The dynamic oscillation threshold y^'j is defined by [5]: 



[ "' \n\d^G[(p{r'-£).r')\dr'^o, 

Jfo+e 



(4) 
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(a) Numerical precision is fixed (prec. = 50). Tj? "^ is computed for 
£=10"* and 3- 10"*. 
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(b) The increase rate of j is fixed (f = 3 ■ 10 *). T^J"" is computed for 
numerical precisions equal to 15 and 100. 



where yo is the initial value of y (i.e. the starting value of 
the linear ramp). Two important remarks can be made 
on this expression (Fig. 6 of [5]): 

• j^J' does not depend on the slope of the ramp e, 
provided that e is small enough, 

• y^'j depends on the initial value yo of the ramp. 

These properties are also observed in numerical simula- 
tions with very high precision. 

2.3 Problem statement 

The above theoretical predictions converge towards the 
observed simulation results for very high numerical pre- 
cision (typically when thousands of digits are considered 



Figure 2: Plot oijiit ^^ a function of the initial condition yo. 
Solid lines are the theoretical prediction y*, ' calculated from 
equation (4). Dashed line represent the values y jj""- 



in the simulation). Figure* 1 shows that for the usual 
double-precision of CPUs (around 15 decimals), theoreti- 
cal predictions of the dynamic bifurcation point y *'' are 
far from thresholds estimated on the numerical simula- 
tion results y j^'"- In particular, the numerical bifurcation 
point y j!"" depends on the slope e, in contrast with the 



*This is a plot similar to figure 10 of [5]. The only difference is that 
the bifurcation point y""™ estimated on the simulation results is here 
defined by the point where the orbit leaves the neighborhood of the 
invariant curve. The motivation for this choice will appear clearly in 
section 3 where random variables are considered. 
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theoretical predictions j^J]. Moreover, figure 2 reveals 
that for a low numerical precision (though even signifi- 
cantly higher than typical precisions used in numerical 
simulations), the dependance of the bifurcation point on 
the initial value 70 is lost over a wide range of 70. 

The minimum precision for which round-off errors do 
not affect the behavior of the system depends on the 
relative magnitude of the slope e, the precision and the 
initial condition Jq. Indeed, figure 1 shows that, beyond 
a certain value of e all curves join the one with highest 
precision. Curves for even higher precisions would over- 
lap, allowing to conclude that they are representative of 
an infinitely precise case. As shown in figure 2, for given 
values of e and of the numerical precision beyond a cer- 
tain value of 70, the theoretical result 7^'j allows to obtain 
a good prediction of the bifurcation delay. 

As a conclusion, the theoretical results obtained in [5] 
are not able to predict the behavior of numerical simula- 
tions carried out at usual numerical precision. The aim of 
this paper is to show how the numerical precision can be 
included in a theoretical model that correctly describes 
numerical simulations. Firstly, it is shown that the model 
computed with a finite precision behaves similarly to the 
model with an ad-hoc additive white noise. This is done 
in the next section. Then, using theoretical results pro- 
posed by Baesens [3], a modified expression describing 
the behavior of the model affected by noise (section 4). 



with an expected value equal to zero (i.e. E[^„] - 0) and a 
level a defined by: 



E[<fm^n] ^CT 5n 



(6) 



where Smn is the Kronecker delta. The definition of the 
expected value E is provided in [14]. For comparison 
with the finite precision case the noise level a is equal to 
10-P'"!. 

The bifurcation point j'l'f"' estimated on the simula- 
tions is defined as the value of 7 for which the orbit leaves 
the neighborhood of the invariant curve. Since the mean 
value of the white noise ^„ is zero, the relevant quantity 
to study is the mean square deviation of the orbit from 
the invariant curve. Therefore, j",'!'" is reached when: 



E[wl]^e, 



(7) 



where Wn- p^- <piYn,£) describes the distance between 
the actual orbit and the invariant curve. Among other 
possible criteria, the condition (7) is chosen because it is 
also used in the analytical calculation made in section 4. 

To simplify the notation, in the rest of the document 
the invariant curve wUl be noted (pij). Its dependency on 
parameter e is no longer explicidy stated. 

In figures 3 and 4, j'^""^ is estimated in the finite preci- 
sion case and in the noisy case. In both cases an average 
is made on w^ obtained in 20 different simulations. Then, 



3. Finite precision versus additive white 

NOISE 

Differences between theoretical predictions and numer- 
ical simulations highlighted in the previous section are 
due to round-off errors that accumulate for finite preci- 
sions. The aim of this section is to verify that round-off 
errors of the computer can be modeled as an additive 
independent and identically distributed random variable 
(referred to as an additive white noise). This result is 
used in section 4 to derive theoretical predictions of the 
dynamic bifurcation point jat- 

Two numerical results are compared. The first is the 
simulation of the system (2) using a numerical precision 
pri (hereafter referred as sl finite precision case). The 
second one (hereafter referred as a noisy case] is the sim- 
ulation of the following stochastic system of difference 
equations: 

Pn^G{p^_i,r„) + in (5a) 

Jn^Jn-i+e, (5b) 

where ^„ is a uniformly distributed stochastic variable 
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Figure 3: Comparison of the dynamic threshold y""™ obtained 
in numerical simulations of a noiseless (2) clarinet model and 
one (5) where a noise of level a = 10~ is introduced. The 
dynamic threshold of oscillation obtained over an average of 
20 runs is plotted against the precision used in the simulations, 
showing that beyond a precision of about a, the system affected 
with noise is insensitive to the precision. 
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Figure 4: Comparison between j'jf" computed for finite precision cases and for noisy cases. For both cases and for each value of e 
we compute the average of the signals Wn = Pn~ 4'iYn) obtained over 20 runs. Then, y"*"" is calculated on the resulting signal. The 
numerical precisions used to simulate the finite precision cases are 7, 15, 30 and 100 decimal digits. 70 = 0. 



Tdl"' ^^ calculated on the mean signal using equation (7). 

In figure 3, J^f" is plotted with respect to the numer- 
ical precision for which systems (2) and (5) are simu- 
lated. The noise level a of the noisy case is equal to 
10"^". Before the numerical precision becomes equal 
to - logjo a - 30, the noise level is smaller than round-off 
errors of the computer. In these situations, the thresholds 
computed in finite precision case and in noisy case are 
equEils. For numerical precisions larger than 30, j'V!'" 
computed on system (5) is constant because the influ- 
ence of the noise overrides that of the round-off errors of 
the computer. Therefore, to avoid the influence of the nu- 
merical precision, the system (5) is now simulated using 
a precision pr2 > 2pri, we choose pr2 - 2pri. 

Figure 4 confirms that the kind of noise introduced in 
the stochastic system can correctly describe a finite preci- 
sion. Indeed, with the exception of the smallest precision 
ipri - 7), the curves are nearly superimposed. Hence, 
in the next section, the stochastic system (5) is studied 
theoretically in order to predict results of numerical sim- 
ulations of system (2) with finite precision. 

4. Analytical study of the noisy dynamic 

CASE 

4.1 General solution of the stochastic clarinet model 

This section introduces a formal solution of the stochastic 
model that is valid when the orbit is close to the invariant 
curve. Function G in equation (5a) is expanded into a 
first-order Taylor series around the invariant curve. Using 
the variable Wn- Pn~ (pijn), the system (5) becomes: 



W„^W„-id^G[(p[j„-£),Yn) + ^n (8a) 

Tn^Tn-i+e. (8b) 

The solution of equation (8a) is [4] : 



w 



„= woY\dxG[(p{jic-e),jt] 



fc=i 



n n 



+ E^fc n dxG[<P^Ym-e],Ym], (9) 

fc=l m=k+l 

where wq is the initial value of w„ . 

Because the additive white noise ^„ has a zero-value 
mean, as in section 3, the relevant indicator is the mean 
square deviation of the orbit from the invariant curve: 
WE[w^]. Equation (9) squared becomes: 



n 



l^\u)Q]\dxG[<p[jk-e),Yk) 



k=\ 

n n 



+ L^fc n dxG[(p[Ym-e),Ym) 

\k=\ m=k+l 



+ 2^0 E n dxG{(p{Yj-£^'Yj)\^k 

n 

X n dxG[(p{Ym-e),jm]. (10) 

m=k+l 

Averaging has no effect on the first term of the right- 
hand side of equation (10) because the stochastic variable 
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Figure 5: Graphical representation of 5xG(p^* (7), 7) and its 
tangent function -1 - ^(7-7$ t) around tlie static oscillation 
threshold for ( = 0.2, 0.5 and 0.8. 



in is not present. Using eq. (6), the average of the second 
term is simplified to: 

'^^L n d^G[ct>ijm-e],jm]\ . (11) 

k=l\m=k+l I 

Because E[(J^„] - 0, the average of the third term of the 
right-hand side of equation (10) is also equal to zero. Us- 
ing the fact that a product can be expressed as an expo- 
nential of a sum of logarithms, the final expression of 
£[14^^] is given by: 



E[w2]^ wl\exp\Y,in\d^G{(pirk-e),jk)\ 



jfc=i 



n 

Y, ln\djcG[(t>ij,n-e),Ym]\ 

m=k+\ 



1\2 



(12) 



Finally, using Euler's approximation, sums are replaced 
by integrals and the expressions of An and B„ become: 

An-wlexpiT" 21n|5;,G((/)(7'-£-),r')|^l, (13) 



1-2 fYn+e 
lya+e 



n~ — 

e Jya 

|exp(|^'^ '21n|a;,G(0(r'-e),r')|^]|dr- 



(14) 



A„ corresponds to the precise case studied in [5] which 
leads to the theoretical estimation j'-J' of the dynamic os- 
cillation threshold for the system without noise (cf. equa- 
tion (4)). Bn is the contribution due to the noise, which 
wiU be considered in the remaining of this section. 

A first glance on equations (13) and (14) allows to ex- 
plain observation made in Section 2.3. Indeed, compar- 
ing the expressions of An and Bn, it possible to distin- 
guish [3, 2] two operating regimes, which, for a given 
values of wq, depends on e, o and 70: 

• An » Bn (deterministic regime) 

In this case the noise does not affect the bifurcation 
delay and the dynamic oscillation threshold can be 
determined by eq. (4). 

• An « Bn (sweep-dominant regime) 

In this case, the bifurcation delay is affected by the 
noise. This regime is studied in the following section. 

In Section 2.3, figures 1 and 2 represent two different 
cases distinguished by the parameter values: in certain 
areas of the figures, the dynamic bifurcation threshold 
does not depend on e but depends on 70, while in other 
areas the dynamic bifurcation threshold depends on e 
but is not dependent on 70. This observation may be 
interpreted as the existence of the two regimes detailed 
above: a sweep-dominant regime and a deterministic 
regime. The transition between the two regimes occurs 
abruptly as observed in figures 1 and 2. 
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Figure 6: Graphical representation ofj'l"'" for different precisions (prec. = 7, 15, 30, 100, 500 and 5000) with respect to the slope e 

and for jo = 0. Results are also compared to analytical static and dynamic thresholds: jst, 7 j! and f'.' (a) j'l'i'^ and ov\y finite 
precision cases are represented, (b) Both finite precision cases and noisy cases are represented for prec. = 7 and 15. 



4.2 Theoretical expression of the dynamic oscillation 
threshold of the stochastic model 

The next step is to find an approximate expression of 



the standard deviation wE [if^] for the sweep -dominant 
regime. In this regime, the term An is negligible with 
respect to the contribution B„ due to the noise, i.e. 
y/e[w^] =; \/B^. The purpose is to obtain a statistical 
prediction of the dynamic oscillation threshold for the 
stochastic system, hereafter referred as f^^. 

It is assumed that e « 1, which implies that the invari- 
ant curve (pij) and the curve p'''*{j] of the fixed points in 
eq. (1) are close [5], and allows the approximation: 



d^G[cl>{r-e),j]^d^G[p^*[r),j]. 



(15) 



Moreover, because of the noise, the bifurcation delay 
is expected to occur earlier, so that the dynamic oscilla- 



tion threshold jiit is assumed to be close^ to the static 
oscillation threshold jst- The term dxG[p^* (71,7) is then 
expanded in a first-order Taylor series around the static 
oscillation threshold 7sf : 



dxG{p^* {j),r) ^ d,G{p^* {Yst),rst) 

=-1 : flip bifurcation 

+ {r-rst)djcyG{p^*{rst),Tst), (i6) 
"^ , ' 

noted -K 

finally we have: 

dxG[p^*{r),r)^-i-K{r-rst), (i?) 

which is used in equation (14). Figure 5 shows the com- 



' This hypothesis could be questioned because according to figures 1 
and 2, even in the presence of noise, the bifurcation delay can be big. 
However, this hypothesis is required to carry out followdng calculations. 
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Using expression (17) the integral 




ni^r' 



/i= [ " '21n|5,G((/)(7'-e),7')| 



contained in the expression (14) of B„ becomes: 



(18) 



h-— [T'-Tst]dr'^-UY-rst] 

E Jv+e t L 



ly+c 



(19) 



The small correction e in the domain of integration can 
be neglected since c « 1. Therefore, we obtain: 



°C^00 0.05 0.10 0.15 0.20 0.25 0.30 

7o 



A = f[(r'-r.Ol 



21 Tn 

r 



(a) Numerical precision is fixed (prec. = 50). y'J and y"™^ are com- 
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(b) The increase rate of 7 is fixed (e = 3- 10"*). f', and y""™2 gj.g com- 
puted for numerical precisions equal to 15 and 100. 



Figure 7: Comparison between theoretical prediction of dy- 
namic oscillation threshold (without noise: j'J] and with noise: 

f*J]] and the dynamic threshold y"""* computed on numerical 
simulations for finite precision case. Variable are plotted with 
respect to the initial condition 70- 



By combining equations (14) and (20), Bn is now writ- 
ten as: 



Bn '■ 



— \ ^■M-\{jn-jstf-{i -jstf\\di 



— exp\-[j„-jst] 



i 



exp — {r'-rst)\dr'. (21) 

7o+e I £ 



The function which appears in the integral I2 is a Gaus- 
sian function with standard deviation 



(22) 



Integral I2 is then [8] 

1 fne 



2Vir''"'lV7(^'-^") 



tTb 



J To 



(23) 



where erf (x) is the error function. The initial condition 70 
is supposed to be much lower than the static threshold 
7st, so that equation (23) can be written: 



parison between dxG[p^* (7), 7) and its tangent function 
-1 - ^^(7- 75f) around the static oscillation threshold for 
( - 0.2, 0.5 and 0.8. The linearisation appears as a good 
approximation of the function in a wide domain of 7 
around the static oscillation threshold jst- For large val- 
ues of the control parameter ( (cf. fig. 5(c)) the linear 
approximation is valid over a narrower range of 7. 



k 



1 fne 



2\ K 



erf a/-(7„-7sO +1 



(24) 



The dependence on the initial condition 70 is now lost. 

Since e « 1, for 7„ > jst the error function quickly 
becomes equal to 1 and the integral I2 is simplified to 
h- v X- Finally the expression of B„ is: 
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Bn '■ 



a 



'M^ir.-r.f 



(25) 



From equation (25) it is possible to obtain the expres- 
sion of \/e [ w^l ~ \/5^: 



which is the theoretical estimation of the dynamic oscil- 
lation threshold of the stochastic systems (5) (or of the 
system (2) computed using a finite precision) when it 
evolves in a sweep-dominant regime. The bifurcation de- 



C,th 



lay is a by-product of eq. (27) since it is simply fj^ - jst 



E[wl]^ae-"'[^]"\xpi^^^{rn-rst) 



(26) 



The dynamic oscillation threshold jHl is defined [3, 2] 
as the value of jn for which the standard deviation 
E [ wf,] leaves the neighborhood of the invariant curve. 



More precisely, the bifurcation occurs when wE [w^] be- 
comes equal to the increase rate e, as defined in eq. (7). 
Finally, using equation (26), we have: 



fdt^rst+\l-^ln 



71^1/1 o 



c-5/4 



(I) 



(27) 



The method presented in this section is based on a 
first-order Taylor expansion of the system (5) around the 
invariant curve </)(7n), leading to the linear system (8). 
Using an asymptotic expansion of the error function it 
is possible to investigate the behavior of \fB^ before y,, 
enters the neighborhood of the static oscillation threshold 
J St- This study allows to define the domain of validity of 
this linear approximation, as done by Baesens [3, 2]. This 
is (T > \fe. This condition is respected in this work since 
a = 10"'"^ with 7<pr< 5000 and 8.10"^ < c < 10"^. More 
details on obtaining the domain of validity are given in 
Appendix A. 
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4.3 Discussion 

In figure 6, f^'j defined by equation (27) is plotted against 
the increase rate e. It is compared with 7 j"*" for different 
values of the precision and for yo = 0. In figure 6(a), 7 jj"" 
is represented for finite precision cases. The significant 
differences between finite precision cases and stochastic 
cases observed for prec. - 7 and 15 are shown in figure 
6(b). The theoretical result f'.'j provides a good estima- 
tion of the dynamic oscillation threshold as long as the 
system remains in the sweep -dominant regime (with a 
better estimation when the bifurcation delay is small*). 
Otherwise, 7^'j is a better approximation of 7^"'", as ex- 
pected in the deterministic regime. 

Figure 7 shows the comparison between f^'j and 7J]"'" 
(only for finite precision cases) plotted against the ini- 
tial condition 70. In figure 7(a), variables are plotted for 
several values of e and for a fixed numerical precision. 
The opposite is done in figure 7(b). As in figure 6, f*,^ 
provides a good estimation of the dynamic oscillation 
threshold in the sweep-dominant regime, as well as 7^^ 
in the deterministic regime. 

Finally, to predict theoretically the dynamic bifurcation 
threshold f^t^ of the stochastic system (5) (as well as of 
the system (2) when it is computed with a finite precision) 
the following procedure is proposed: 

• compute the theoretical estimation f^'j of the 
stochastic system through eq. (27) 

• compute the theoretical estimation 7^^ of the system 
without noise through eq. (4) 

• if 7^'j < 7^'j the system remains in the "sweep - 
dominant regime" and the dynamic threshold f^'j is 
equal to 7'?*, otherwise the "deterministic regime" is 
attained and the dynamic threshold 7^^ is equal to 

th 
' df 

Figure 8 compares the relative error RE of the three 
theoretical predictions of the oscillation threshold (7^^, 
7'/' and f'/p with respect to 7j!"". as a percent value: 



RE{X\ = 100 X 



i7gr-^i 

' dt 



(28) 



where X takes successively the values of 7st, 7^^^ and 7^^. 
For standard double-precision (fig. 8(a), prec.=15), the 
sweep-dominant regime is prevalent throughout most of 
the range of increase-rates studied in this article. Higher 
precisions (for instance prec.=100) imply the appearence 



'^This is an expected result because of the initial assumption of a 
small bifurcation delay in the presence of noise, leading to first-order 
Taylor expansions jst in previous calculation (see step one in section 
4.2). 



of the deterministic-regime for lower increase-rates. In 
this case, fl'j provides a better estimation of the oscilla- 
tion threshold of the clarinet with a linearly increasing 
blowing pressure. Indeed, in situations represented in 
figure 8, RE [7 j'!] never exceeds 15% while RE [jst] and 
RE [jlli] can reach 60% and 145% respectively. At slightly 
lower values of e than the limit between the two regimes, 
j'J] still provides a better estimation of j'i'f'" than f^Jl, a 
situation that occurs for all values of the precision, ac- 
cording to figure 6. 

5. Conclusion 

In many situations, the finite precision used in numerical 
simulations of the clarinet system does not produce ma- 
jor errors in the final results that are sought. Such is the 
case, for instance, when estimating the amplitudes for a 
given regime. 

However, when slowly increasing one of the control 
parameters, the distances between the state of the system 
and the invariant curve can become smaller than the 
round-off errors of the calculation, with dramatic effects 
on the time required to trigger an oscillation. In these 
cases, the inclusion of a stochastic variable in the theory 
allows to correctly estimate the threshold observed in 
simulations, which lies between the static and dynamic 
thresholds found for precise cases. 

As a final remark, the present theoretical study is prob- 
ably not restricted to describe numerical simulations. In- 
deed, the noise level a measured in an artificially blown 
instrument is typically of the order of magnitude of 10"^. 
The domain of validity of the results: c > \/e suggests 
that the comparison with experiment using blowing pres- 
sure with increase rates e > 10~^, seems to be possible 
although the noise level usually increases with the pres- 
sure applied to the instrument. 
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A. Limit of the linear calculation 

The method presented in Section 4 is based on a first- 
order Taylor expansion of the system (5) around the in- 
variant curve (pijn) leading to define the linear system (8). 
Following Baesens [3, 2], we give here the upper bound 
of the domain of validity of this linear approximation. 
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Using equations (14) and (24), the expression B„ is 
given by: 



Br, — 



;,»p|^,r„-.,)- 



erf h/-(r«-rsf) +^ 



(29) 



We investigate the behavior of E [ if ^] before 7„ enters 
in the neighborhood of the static oscillation threshold 
jst- More precisely, we compute an approximate expres- 
sion of E [ If ^j] when jn < Yst ~ 'v, where v is defined by 
equation (22) . To do this, the error function in equation 
(29) is expanded in a first-order asymptotic series [1] (the 
asymptotic expansion of the error function erf (x) for large 
negative x is recalled in Appendix B): 



0-2 nr IK, ,2 



exp(-f (7„-7,f)^] 
-1 ^^ '- + 1 



[jn-Tst) 



(30) 



which is simplified in: 



(31) 



2K[rn-rstY 

Using the explicit form of 7„, solution of equation (2b): 



jn^en + jQ, 



and (31), we have: 



'Bn - 



1 



(32) 



(33) 



^/2Ke\/nst-n 
where Ust is the iteration for which jst is reached. 

Equation (33) means that when Yn<Tst- f^, the stan- 
dard deviation WE [if^] =; \/B^ increases with the time 
(i.e. with n) like \l ^Ust- n to order a/^, and there- 
fore remains small if a « \/c. Otherwise, if cr > v^, the 
orbit of the series p^ leaves the neighborhood of invari- 
ant curve <p(j) before the static oscillation threshold is 
reached. In this case, linear calculation made in Section 
4 is no longer valid. 

B. Asymptotic expansion of error 

FUNCTION 

The asymptotic expansion of the error function erf(x) for 
large negative x (x ^ -oo) is [1]: 



erf(x) : 



exp(-x^) 

x/ttx 



1+ E (-1)' 

m=l 



1-3... (2m -1) 

(2x2)'" 



(34) 
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